masterData <- read_csv("analysis_data.csv")

set.seed(123)

library(stargazer);library(sjPlot);library(gridExtra);library(ggrepel)
library(arm); library(scales);library(lme4)




#TABLE 4: Non-Nested Multi-Level Models with Random Intercepts: Determinants of Liberal Votes




mIdeology<-glmer(direction ~    scIdeolScore +
                       (1|ChiefJustice)+ (1|HCDBcaseId)+
                       (1|primaryIssueSubArea)+(1|justice),
                     data=masterData[masterData$unan==0,], family = 'binomial',
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

mControlsIdeology<-glmer(direction ~    scIdeolScore +
                           practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                           practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                   (1|ChiefJustice)+ (1|HCDBcaseId)+
                   (1|primaryIssueSubArea)+(1|justice),
                 data=masterData[masterData$unan==0,], family = 'binomial',
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

mControlsIdeologyLC<-glmer(direction ~ scIdeolScore+ lcDispositionDirection +lcDissent+
                        practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                          practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                        (1|ChiefJustice)+ (1|HCDBcaseId)+
                        (1|primaryIssueSubArea)+(1|justice),
                      data=masterData[masterData$unan==0,], family = 'binomial',
                      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))




#Regression table
stargazer(mIdeology, mControlsIdeology, mControlsIdeologyLC,
          type="html", align=T,style="ajps",
          star.cutoffs = c(.05, .01, .001),
          out='Manuscript Regression Table 1.html')





#FIGURE 1: Predicted Probabilities With Judge-Level Ideology Score, All Non-Unanimous Cases 


#Creating a sample data set with a row for each justice and covariates based on the mean values for each judge.



#getting the model dataset (same observations for all models; 
#this frame includes appPMParty, which is used in the visualisation)


model_data<-mControlsIdeologyLC@frame



# a function to get the modal value of any vector
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}


# adding information to inform random effects prediction
# includes justice; primaryIssueSubArea;Chief Justice
# It excludes the ID of case for parsimony

model_data$justice<- as.character(model_data$justice)

model_data<- model_data%>% sjmisc::to_dummy(justice, var.name = "V", suffix = "numeric") %>% 
  bind_cols(model_data) 

model_data<- model_data%>% sjmisc::to_dummy(primaryIssueSubArea, suffix = "numeric") %>% 
  bind_cols(model_data) 

model_data<- model_data%>% sjmisc::to_dummy(ChiefJustice, suffix = "label") %>% 
  bind_cols(model_data) 

model_sum<-model_data%>%group_by(justice)%>%
  summarise(scIdeolScore=getmode(scIdeolScore),
            lcDispositionDirection   = mean(as.numeric(lcDispositionDirection=='liberal')),
            lcDissent = mean(lcDissent),
            practiceSpecEconomic = getmode(practiceSpecEconomic), 
            practiceSpecCriminal = getmode(practiceSpecCriminal), 
            practiceSpecPublic = getmode(practiceSpecPublic), 
            practiceSpecCivil = getmode(practiceSpecCivil), 
            practiceSpecCommon = getmode(practiceSpecCommon), 
            female = getmode(female), 
            PriorExp = getmode(PriorExp),
            ChiefJustice= getmode(ChiefJustice),
            HCDBcaseId= getmode(HCDBcaseId),
            primaryIssueSubArea= getmode(primaryIssueSubArea),
            meanDirection= mean(direction),
            justice33=mean(justice_1),
            justice34=mean(justice_2),
            justice35=mean(justice_3),
            justice36=mean(justice_4),
            justice37=mean(justice_5),
            justice38=mean(justice_6),
            justice39=mean(justice_7),
            justice40=mean(justice_8),
            justice41=mean(justice_9),
            justice42=mean(justice_10),
            justice43=mean(justice_11),
            justice44=mean(justice_12),
            justice45=mean(justice_13),
            justice46=mean(justice_14),
            justice47=mean(justice_15),
            justice48=mean(justice_16),
            justice49=mean(justice_17),
            justice50=mean(justice_18),
            justice51=mean(justice_19),
            justice52=mean(justice_20),
            justice53=mean(justice_21),
            subArea1=mean(primaryIssueSubArea_1),
            subArea2=mean(primaryIssueSubArea_2),
            subArea3=mean(primaryIssueSubArea_3),
            subArea4=mean(primaryIssueSubArea_4),
            subArea5=mean(primaryIssueSubArea_5),
            subArea6=mean(primaryIssueSubArea_6),
            subArea7=mean(primaryIssueSubArea_7),
            subArea8=mean(primaryIssueSubArea_8),
            subArea9=mean(primaryIssueSubArea_9),
            subArea10=mean(primaryIssueSubArea_10),
            subArea11=mean(primaryIssueSubArea_11),
            subArea12=mean(primaryIssueSubArea_12),
            subArea13=mean(primaryIssueSubArea_13),
            subArea14=mean(primaryIssueSubArea_14),
            subArea15=mean(primaryIssueSubArea_15),
            subArea16=mean(primaryIssueSubArea_16),
            subArea17=mean(primaryIssueSubArea_17),
            subArea18=mean(primaryIssueSubArea_18),
            subArea19=mean(primaryIssueSubArea_19),
            subArea20=mean(primaryIssueSubArea_20),
            subArea21=mean(primaryIssueSubArea_21),
            subArea22=mean(primaryIssueSubArea_22),
            subArea23=mean(primaryIssueSubArea_23),
            subArea24=mean(primaryIssueSubArea_24),
            subArea25=mean(primaryIssueSubArea_25),
            subArea26=mean(primaryIssueSubArea_26),
            subArea27=mean(primaryIssueSubArea_27),
            subArea28=mean(primaryIssueSubArea_28),
            subArea29=mean(primaryIssueSubArea_29),
            subArea30=mean(primaryIssueSubArea_30),
            subArea31=mean(primaryIssueSubArea_31),
            subArea32=mean(primaryIssueSubArea_32),
            subArea33=mean(primaryIssueSubArea_33),
            Brennan=mean(ChiefJustice_Brennan),
            French=mean(ChiefJustice_French),
            Gleeson=mean(ChiefJustice_Gleeson),
            Keifel=mean(ChiefJustice_Keifel),
            numberCases= n())







model_sum$initials<-judge_data$justiceInitials[1:21]
model_sum$appPMParty<-judge_data$appPM[1:21]

model_sum$appPMParty[model_sum$appPMParty==2]<- "coalition"
model_sum$appPMParty[model_sum$appPMParty==1]<- "alp"

# predicted probability function for the ideology model
mfxIdeologyModel<- function(model,newdata){
  sim <- sim(model, 1000)
  
  newdata$pr_mean<-NA
  newdata$pr_lower<-NA
  newdata$pr_upper<-NA
  newdata$pr_sd<-NA
  
  for(i in 1:nrow(newdata)){
    #get the  point estimate
    predictJustice<- sim@fixef[,1]+
      sim@fixef[,2]*newdata$scIdeolScore[i]+
      sim@fixef[,3]*newdata$lcDispositionDirection[i]+
      sim@fixef[,4]*newdata$lcDissent[i]+
      sim@fixef[,5]*newdata$practiceSpecEconomic[i] +
      sim@fixef[,6]*newdata$practiceSpecCriminal[i] +
      sim@fixef[,7]*newdata$practiceSpecPublic[i] +
      sim@fixef[,8]*newdata$practiceSpecCivil[i] +
      sim@fixef[,9]*newdata$practiceSpecCommon[i] +
      sim@fixef[,10]*newdata$PriorExp[i] +
      sim@fixef[,11]*newdata$female[i]


    
    
    #confidence intervals
    newdata$pr_mean[i] <- invlogit(mean(predictJustice))
    newdata$pr_lower[i]   <- quantile(invlogit(predictJustice), 0.025)
    newdata$pr_upper[i] <- quantile(invlogit(predictJustice), 0.975)
    newdata$pr_sd[i] <- sd(invlogit(predictJustice))  
    
  }
  return(newdata)
  
}



#applying predicted values function
model_sumIdeology<-mfxIdeologyModel(mControlsIdeologyLC,model_sum)

#plot

plot_ideol<-ggplot(model_sumIdeology, aes(scIdeolScore, pr_mean))+ 
  geom_point()+
  geom_linerange(aes(ymin=pr_lower, ymax=pr_upper))+ 
  scale_y_continuous(limits = c(0,1))+geom_smooth(method = 'lm', se=F, col='Black')+
  scale_x_continuous(limits = c(0,1))+
  geom_text_repel(aes(scIdeolScore, pr_mean, label=rev(model_sum$initials), col=rev(model_sum$appPMParty)),nudge_x=0.01, show.legend = F)+
  theme_classic()+xlab("Liberal ideology")+ylab('pr(Liberal vote)')+
  scale_color_manual(values=c('red', 'royalblue3'))

plot_ideol




ggsave(plot_ideol,filename = 'figure_1_Colour.png' ,height = 5, width = 7, device='png', dpi=600)





## figure 1 (black and white)


plot_ideol<-ggplot(model_sumIdeology, aes(scIdeolScore, pr_mean))+ 
  geom_point()+
  geom_linerange(aes(ymin=pr_lower, ymax=pr_upper))+ 
  scale_y_continuous(limits = c(0,1))+geom_smooth(method = 'lm', se=F, col='Black')+
  scale_x_continuous(limits = c(0,1))+
  geom_text_repel(aes(scIdeolScore, pr_mean, label=rev(model_sum$initials), col=rev(model_sum$appPMParty)),nudge_x=0.01, show.legend = F)+
  theme_classic()+xlab("Liberal ideology")+ylab('pr(Liberal vote)')+
  scale_color_manual(values=c('Black', 'grey'))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA))

plot_ideol



ggsave(plot_ideol,filename = 'figure_1_BW.png' ,height = 5, width = 7, device='png', dpi=600)




#### Interaction model.

mControlsIdeologySpecPracticeVaryingSlopes <-glmer(direction ~  scIdeolScore*Area+
                                                     practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                                                     practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                                                     lcDispositionDirection +lcDissent+
                                                     (1|ChiefJustice)+ (1|HCDBcaseId)+
                                                     (1|primaryIssueSubArea)+(1|justice),
                                                   data=masterData[masterData$unan==0 & masterData$Area!="International and Maritime Law",], family = 'binomial',
                                                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
summary(mControlsIdeologySpecPracticeVaryingSlopes)


# manuscript table 5
tab_model(mControlsIdeologySpecPracticeVaryingSlopes, show.est= T,show.re.var	=T,transform = NULL, show.r2	=F, 
          file= "Manuscript Regression Table 2.html")

plot_model(mControlsIdeologySpecPracticeVaryingSlopes, type = "int" ,terms='scIdeolScore', grid=T,
           colors = "black")+
  ggplot2::scale_y_continuous(limits = c(0,1))+ 
  ggplot2::scale_x_continuous(limits = c(0,1))+ggtitle("")+
  ggplot2::xlab("Liberal Ideology")+ggplot2::ylab('pr(Liberal vote)')+theme_classic()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA))
ggplot2::ggsave(filename = 'figure_2.png' ,height = 6, width = 8, device='png', dpi=600)


#manuscript table 7

library(brms)



mod <- brm(direction ~ me(scIdeolScore, sdBoot) + 
             practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
             practiceSpecCivil+practiceSpecCommon+PriorExp+female+
             lcDispositionDirection + lcDissent + 
             (1|HCDBcaseId) + (1|primaryIssueSubArea)+ 
             (1|ChiefJustice)+ (1|justice),
           data=model_data,
           family = bernoulli)
summary(mod)